#############################################
#############################################
### REPLICATION MATERIAL                  ###  
### Gender Quotas and Political Knowledge ###
### Women in Parliament as main IV (main) ###
#############################################
#############################################

# Libraries needed
library(tidyverse)
library(stargazer)
library(haven)
library(lme4)
library(rstanarm)
library(tidybayes)
library(BayesPostEst)
library(ggplot2)
library(ggridges)
library(magrittr)
library(dplyr)
library(purrr)
library(forcats)
library(tidyr)
library(modelr)
library(ggdist)
library(tidybayes)
library(cowplot)
library(rstan)
library(RColorBrewer)
library(gganimate)
library(wrangle)

# Set plot theme
theme_set(theme_tidybayes() + panel_border())

# Set seed and stan options (to make it run faster)
set.seed(123)
options(mc.cores = parallel::detectCores())

# Read RDS data
analysis.2 <- readRDS("data/analysis/analysis_2_df.rds")

# Modify full dataset
analysis.int <- analysis.2 %>%
  mutate(female = recode(gender,
                         "1" = 0,
                         "0" = 1)) %>%
  select(ID, country, name, year, know_q_num, know_correct, women_in_parl,
         female, quota, age, education, GDP_growth, counter) %>%
  mutate(age = cut(age, 5,labels = FALSE)) %>%
  mutate(education = cut(education, 5, labels = FALSE)) %>%
  mutate(treated = ifelse(name == "Belgium" | name == "France" | name == "Greece" |
                            name == "Ireland" | name == "Portugal" |
                            name == "Spain" , 1, 0)) %>%
  write_rds("data/analysis/analysis_int.rds") %>%
  glimpse()


# Binomial version of the dataset - full sample
analysis.int.binomial <- analysis.int %>%
  group_by(female, quota, age, treated, counter, education, women_in_parl,
           GDP_growth, country, year, know_q_num) %>%
  summarize(n_correct = sum(know_correct),
            n_incorrect = n() - sum(know_correct)) %>%
  write_rds("data/analysis/analysis_int_binomial.rds") %>%
  glimpse()

###################################
###################################
## MODEL AND RESULTS TABLE (A4) ###
###################################
###################################

# Model
fit.women <- stan_glmer(cbind(n_correct, n_incorrect) ~ female*age*women_in_parl +
                          education + 
                          GDP_growth + 
                          (1 | country) + 
                          (1 | year) + 
                          (1 | know_q_num) + 
                          (1 + female | country:know_q_num), 
                        data = analysis.int.binomial,
                        family=binomial(link="logit"),
                        iter = 6000,
                        control = list(adapt_delta = 0.99))

#######################
### Summary - TABLE A4
#######################
print(fit.women, digits = 4)
mcmcReg(fit.women, format = "latex")
mcmcTab(fit.women)

#######################
### Convergence checks
#######################

# Rhat - convergence ok
dev.copy(png,'rhat_fit_int.png')
plot(fit.women, "rhat")
dev.off()

# Prior summary
prior_summary(fit.women)

#################################
### PREDICTED PROBS DATASETS  ###
#################################

# Dataframe with predicted probabilities (no random effects)

pred_fit.women_noRE <- crossing(country = unique(analysis.int.binomial$country), 
                                women_in_parl = 4.7:41.10, 
                                female = c(0,1),
                                know_q_num = c(1:10),
                                age = c(1:5),
                                quota = c(0,1),
                                education = mean(analysis.int.binomial$education, 
                                                 na.rm = T),
                                GDP_growth = mean(analysis.int.binomial$GDP_growth, 
                                                  na.rm = T)) %>%
  add_epred_draws(fit.women, 
                  scale = "response", 
                  re_formula = NA) %>%
  write_rds("data/analysis/pred_fit_women_noRE.rds")

# Dataframe with draws
draws <- fit.women %>%
  spread_draws(b[term,group]) %>%
  filter(str_detect(group, "country:know_q_num:")) %>%
  separate(group, c("country_2", "q_num_2", "country", "q_num"), ":") %>%
  write_rds("data/analysis/draws_fit_women.rds")

# Dataframe with predicted probabilities (with random effects)
pred_fit.women_c_RE <- crossing(country = unique(analysis.int.binomial$country), 
                                women_in_parl = 4.7:41.10, 
                                female = c(0,1),
                                year = 2005,
                                know_q_num = c(1:10),
                                age = c(1:5),
                                education = mean(analysis.int.binomial$education, 
                                                 na.rm = T),
                                GDP_growth = mean(analysis.int.binomial$GDP_growth, 
                                                  na.rm = T)) %>%
  add_epred_draws(fit.women, 
                  scale = "response", 
                  re_formula = ~ (1 + female | country:know_q_num)) %>% 
  write_rds("data/analysis/pred_fit_women_c_RE.rds")

#############
#############
### PLOTS ###
#############
#############


##################################################################
### Figure 6: Predicted Probabilities - Women Parl and Gender  ###
##################################################################  

plot_women_1 <- readRDS("data/analysis/pred_fit_women_noRE.rds")

plot_women_1$female <- factor(plot_women_1$female,
                              labels = c("Male", "Female"))
plot_women_1$quota <- factor(plot_women_1$quota, 
                             labels = c("No Quotas", "Quotas"))
plot_women_1$age <- factor(plot_women_1$age, 
                           labels = c("15-31", "32-48", "49-65", "66-82", "83-99"))

plot_women_1 %>%
  ggplot(aes(x = women_in_parl, y = .epred, fill = female, color = female)) +
  stat_pointinterval() +
  facet_grid(female ~ .) +
  scale_color_grey(name = "Gender",
                     labels = c("Male", "Female")) +
  scale_fill_grey(name = "Gender",
                    labels = c("Male", "Female")) +
  labs(x = "Women in Parliament", y = "Predicted Probabilities")
ggsave("plots/main/fig-6-greyscale.png")

#####################################################################
### Figure 6: Predicted Probabilites - Women Parl and Gender/Age  ###
##################################################################### 

plot_women_1 %>%
  ggplot(aes(x = women_in_parl, y = .epred, color = female, fill = female)) +
  stat_pointinterval() +
  facet_grid(factor(age, levels=rev(unique(age))) ~ female) +
  scale_fill_grey(name = "Gender",
                     labels = c("Male", "Female")) +
  scale_color_grey(name = "Gender",
                     labels = c("Male", "Female")) +
  labs(x = "Women in Parliament", y = "Predicted Probabilities")
ggsave("plots/main/fig-7-greyscale.png")

#####################################################################
### Figure A3: Percentage of Women in Parliament by Country       ###
##################################################################### 

women_parl <- analysis.2 %>%
  group_by(year) %>%
  select(name, year, women_in_parl, quota) %>%
  ungroup() %>%
  distinct(name,year,women_in_parl,quota) %>%
  write_rds("data/appendix/women-parl.rds") %>%
  glimpse()

women_parl %>%
  ggplot(aes(x = year, y = women_in_parl)) +
  geom_line() + 
  facet_wrap(. ~ name) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) +
  theme_bw() + 
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()) +
  labs(x = "Year",
       y = "Percentage of Women in Parliament") +
  geom_vline(data=filter(women_parl, name == "Belgium"), 
             aes(xintercept = 1999), colour="red") +
  geom_vline(data=filter(women_parl, name == "France"), 
             aes(xintercept = 2000), colour="red") +
  geom_vline(data=filter(women_parl, name == "Greece"), 
             aes(xintercept = 2012), colour="red") +
  geom_vline(data=filter(women_parl, name == "Portugal"), 
             aes(xintercept = 2006), colour="red") +
  geom_vline(data=filter(women_parl, name == "Spain"), 
             aes(xintercept = 2007), colour="red") +
  geom_vline(data=filter(women_parl, name == "Ireland"), 
             aes(xintercept = 2012), colour="red")

ggsave("plots/appendix/women-parl.png")



